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Abstract. We introduce a model for a population on a lattice with diffusion and 
birth/death according to 2 A — ► 3A and A — ► <fi for a particle A. We find that 
the model displays a phase transition from an active to an absorbing state which is 
continuous in 1 + 1 dimensions and of first-order in higher dimensions in agreement 
with the mean field equation. For the 1 + I dimensional case, we examine the critical 
exponents and a scaling function for the survival probability and show that it belongs 
to the universality class of directed percolation. In higher dimensions, we look at the 
first-order phase transition by plotting a histogram of the population density and use 
the presence of phase coexistence to find an accurate value for the critical point in 
2 + 1 dimensions. 



PACS numbers: 64.60.Ht, 05.70.Jk 
1. Introduction 

Non-equilibrium phase transitions have long been a major area of investigation (see [1] 
for a review). Studies are wide-ranging, involving atmospheric precipitation [2], sandpile 
models [3], epidemics [I] and many more. One of the great achievements in the field 
has been the discovery that a broad range of models belong to one of a few universality 
classes, whose members share the same critical exponents and scaling functions (see 
for example [5]). In particular, it has been conjectured that all models with a scalar 
order-parameter that exhibit a continuous phase transition from an active state to a 
single absorbing state belong to the same universality class of directed percolation (DP) 
[HI Ej. Although not proven, the conjecture is strongly supported by numerical evidence 
and seems to be even more general since, for example, a system with multiple absorbing 
states is known to belong to the class [8]. 

Here, we introduce a lattice model representing a population in a habitat and 
include the processes of birth, death and diffusion. Due to the conflict between growth 
and decay, with steady state population density p as our order parameter, we expect a 
phase transition to an absorbing state to occur under certain conditions. The expected 
phase transition may be either continuous or of first-order. The latter, however, is rarely 
seen in low spatial dimensions due to the destabilisation of the ordered phase caused by 
the larger fluctuations that are present in such systems. Hinrichsen has hypothesised 
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that first-order phase transitions are impossible in 1 + 1 dimensional systems provided 
that there are no additional conservation laws, long-range interactions, macroscopic 
currents or special boundary conditions [9]. Since our model does not fulfil any of these 
criteria, we expect at least the 1 + 1 dimensional version of our model to exhibit a 
continuous phase transition and thus to belong to DP. 

Thanks to series expansions, the critical exponents of DP are now known to a high 
degree of accuracy [10]. We therefore proceed in Section [2] by describing our model 
and then, in Sections [3] and HI examine the critical exponents and a scaling function 
respectively to compare them with those of DP. In Section [5l we examine the first-order 
phase transitions and conclude with some remarks in Section El 



2. The Model 



We have a d- dimensional square lattice of linear length L where each square is either 
occupied by a single particle (1) or is empty (0). A site is chosen at random. With 
probability pd the particle on an occupied site dies, leaving the site empty. If the particle 
does not die, a nearest neighbour site is randomly chosen. If the neighbouring site is 
empty the particle moves there; otherwise, the particle reproduces with probability p^ 
producing a new particle on another randomly selected neighbouring site, conditional 
on that square being empty. A time step is defined as the number of lattice sites 
N = L d and periodic boundary conditions are used. We have the following reactions 
for a particle A for proliferation and annihilation respectively, 

A + A — >3A and A — (1) 

Our model is similar to that of Schlogl's second model pj] except for the inclusion of 
diffusion and the reactions given in ([1]) being unidirectional. 

Assuming the particles are spaced homogeneously, the mean field equation for the 
density of active sites p(t) is given by 

^ = Pb(l - Pd)p(t) 2 (l - p{t)) - p dP (t). (2) 
This has three stationary states 



Po = 0, p± = ~(l±«/l- 4pd ). (3) 
2 V V Pb(l - Pd)/ 

Clearly, for 4p d > p h (l — p d ), p is the only real stationary state, resulting in a phase 

transition occurring at the critical death rate p^ c = Pb/(4 + Pb)- Simple analysis 

shows that p + and po are stable stationary states, whereas p_ is unstable and therefore 

represents a critical density p c below which extinction will occur in all cases. So, for 

Pd < Pd c , 

/,\ f for pit) < Pc , /,\ 

p(t) — ► <^ _ '){ yc as t — ► oo. 4 

I P+ for p(t) > p c 

At pa = Pa c1 the stationary density jumps from 1/2 to 0, resulting in a first-order 
phase transition. We investigate whether the Monte Carlo (MC) simulations agree with 
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the mean field by plotting, in Figure [fl the steady state population density against p d 
in 1 + 1, 2 + 1 and 3 + 1 dimensions with p(0) = 1. From now on, we use a constant 
value of pb = 0.5. We note that, although we keep p\, constant, the actual birth rate is 
dependent on p d since the probability of birth is proportional to pt>(l — Pd)- To find the 
steady state, we examine surviving runs only, looking at an increasing number of time 
steps up to 2 x 10 5 and increments in pd of 5 x 10~ 5 as the critical point is approached. 

We see from the results that we have a strong indication of a first order phase 
transition in 2 + 1 and 3 + 1 dimensions whilst a continuous phase transition in 1 + 1 
dimensions. While these results are compelling, we note that since the simulations were 
performed on finite lattices and for finite times, we cannot take them to be conclusive 
since in such simulations there is always a non-zero probability of survival for finite t 
even for pd > Pd c - Instead we look for power-law behaviour in p(t) close to the critical 
point. For a continuous phase transition we expect asymptotic power-law behaviour of 
the order parameter at the critical point (see [12] for a review) of the form 

p{t) oc r s . (5) 

In log-log plots, positive curvature for large t indicates the system is in the active phase 
whereas negative curvature implies that the system is in the absorbing phase. A first 
order phase transition will therefore be marked by non power-law behaviour, rather 
exponential decay of the order parameter for pd > Pd c - Figure [2] shows the 1 + 1 
dimensional case with p(t) for different values of pd close to the critical point clearly 
showing power-law behaviour at the critical point. No power-law behaviour is however 
observed in 2 + 1 or 3 + 1 dimensions; instead, with the inclusion of spontaneous particle 
creation at rate k = 0.005, hysteresis occurs in both cases as plotted in Figure [31 which 
is indicative of first-order phase transitions. However, no hysteresis is observed in the 
1 + 1 dimensional case. We examine the first-order phase transitions in more detail in 
Section [51 
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Figure 2. Power-law behaviour for the 1 + 1 dimensional model. Solid lines represent 
(from top to bottom) p d = 0.071654, 0.071754 and 0.071854. The hashed line 
represents the gradient -0.159 as a guide for the eye. The insert shows non-power law 
behaviour for various values of pd close to the critical point for the 2 + 1 dimensional 
model. The 3 + 1 dimensional case is very similar. Information on how the critical 
points were found are detailed later in the paper. 
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Figure 3. Hysteresis loop for the (a) 2 + 1 and (b) 3+1 dimensional models with p& 
increasing (x) and decreasing (+). The insert in (a) shows no hysteresis occurring in 
the 1 + 1 dimensional model. 



Since the model shows a continuous phase transition in 1 + 1 dimensions, we must 
ask ourselves to which universality class it belongs. By Grassberger and Janssen's 
conjecture [6j [7], we would expect it to belong to the universality class of DP. We 
investigate this now by looking at the critical exponents and a scaling function in turn 
for 1 + 1 dimensions only. 
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3. Critical Exponents 

Finding the critical exponents through steady state simulations is notoriously difficult 
due to critical slowing down, finite-size effects, large fluctuations and the difficulties 
that arise in finding the critical point. A much more effective method is that of time- 
dependent simulations, which has proved to be a very efficient way of determining the 
critical exponents and the critical point for models exhibiting absorbing phase transitions 
[13] . Using this method, the time-evolution of the model is observed up to some time 
£mj after beginning with a configuration that is very close to the absorbing state - 2 
adjacent particles in this model. The size of the lattice is made large enough so that 
the particles never reach the boundary before tu- 

We measure the survival probability P(t), defined as the probability that the system 
has not reached the absorbing state at time t and the average number of occupied sites 
n{t). At pd = pd C ) we expect the following asymptotic power-law behaviour [T4] 

P(t) oc t' 6 , (6) 
n(t) cxf. (7) 

Away from the critical point, the evolution departs from pure power-law and so by 
examining log-log plots of P(t) and n(t) versus t we can find the critical point by finding 
the value of pd that gives a straight line. Generally, however, we expect corrections to 
the pure power-law behaviour so that P(t) is more accurately given as [T4"j 

P(t) ocr s (l + ar 1 + br s ' + (8) 

and similarly for n(t). Here, 5' represents the correction-to-scaling exponent for 5. More 
precise estimates for the critical exponents are obtained by examining the local slope 

_ 

m(m) 

and similarly for n(t), where the critical exponent 5 is given by lim^oo 5(t). Here, m is 
given as the local range over which the slope is measured and is typically 5 [T5l [16] or 
8 [H]. Grassberger [H] has shown that for the local slope defined in we have the 
following behaviour 

5(t) = 5 + at- 1 + b5't- & ' + ... (10) 

and again, similarly for n(t). Thus if we plot 5(t) versus we have that the critical 
exponent 5 is given by the intercept with the y axis and any curvature would indicate 
a correction-to-scaling exponent less than 1. 

We plot in Figure HJ v{t) an d S(t) for up to t = 10 6 and over 10 5 runs. We find that 
the data is in fact very noisy, leading to inaccurate results for the critical exponents. 
However, it is clear that the gradient increases for pd = 0.071746 and decreases for pd = 
0.071762 for large t whereas we have an approximately straight line for pd = 0.071754, 
especially for rj(t). This therefore gives a value of pd c = 0.071754 ±0.000004. We plot in 
Figure [5] i](t) and 5(t) for p d = 0.071754 but this time up to t = 10 4 and over 2.5 x 10 6 
runs to improve the accuracy of the results. From this plot we can read off the values 
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Figure 4. Plots of (a) r)(i) and (b) — 5(t) up to t 
p d = 0.071746, 0.071754 and 0.071762. 
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Figure 5. Plots of (a) rj(t) and (b) -6(t) up to t = 10 4 at p d = p dc = 0.071754. 



7] = 0.312 ±0.002 and 5 = 0.160 ±0.002, which are in agreement with the best currently 
known DP values of 77 = 0.313686 and 5 = 0.159464 [10]. 



4. Scaling Functions 

For models to belong to the same universality class, both the critical exponents and the 
scaling functions must be the same. Having shown that the critical exponents S and rj 
for our model are in agreement with the DP values, we turn our attention now to the 
scaling function for the probability of survival, P(A, L, t), where A = p^ c — pd- We have 
the scaling ansatz [5] 

P{A,L,t) ~ \- f> R 1 (A\,L\-'' J ;t\- u i) (11) 
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for any A > 0. (3 is the critical exponent associated with the survival probability 
according to 

P ~ A? (12) 

as t — > oo and un and u± are the critical exponents associated with the temporal 
and spatial correlation lengths respectively. As in the previous section, we begin our 
simulations from a single seed with L large enough so that we can ignore the in- 
dependence in ffTTl) and examine 

P(A, t) ~ \- p R 2 {A\, t\-"u), (13) 

only. The functions Ri are non-universal, i.e. they are unique to each model. However, 
we may introduce metric factors a a and a t , 

P(A, *) ~ \-PR(a A A\, otU - "" ) (14) 

so that R is a universal scaling function. Now all of non- universal, system-dependent 
features, such as the update scheme, boundary conditions, lattice structure, etc. are 
contained in these non-universal metric factors. 

If we choose A = (a^) 1 '^, then, since 5 = f3/v\\, 

P(A, t) ~ (a t t)- <5 P(a A A(a t t) 1/iy n , 1) (15) 

so we would expect to observe a data collapse by plotting (a t t) 5 P versus a A A(aib) l ' v w 
for all models belonging to the same universality class. To find the metric factors a t 
and a a, we use the normalisations 

P(l,oo) = P(0, 1) = 1. (16) 
Then, choosing A = and A = {a t t) 1 ^^ we have that a t is given by the amplitude of 

P~(a t ty s . (17) 
Similarly, setting A = (a^A)' 1 and letting t — ► oo, a A is given by the amplitude of 

P ~ (a A Af, (18) 

where, since A > 0, we must have A > 0. 

We now compare the scaling function for our model with that for directed bond 
percolation which belongs to DP. The model begins with an initial number of active sites 
and proceeds in time with site i becoming active at time t + 1 if site i + 1 and/or site 
z — 1 is active at time t and there exists a bond (with probability p) between this active 
site and the site i. The best estimate of the critical value of p is p c = 0.644700185(5) 
[T7] where, for p > p c , the active sites percolate the infinite system. 

As was the case in our model, we begin with a single seed - just one active site at 
the origin here - and record the probability of survival Pop(t), i.e. the probability that 
there is at least one active site at time t. We denote the metric factors with a subscript 
DP for the directed bond percolation model and subscript P for our model. Excluding 
the metric factors, we plot in the insert in Figure [6] the two separate data collapses for P 
for both the directed bond percolation model and for our own by plotting Pt s vs. tA"n . 
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Figure 6. Data collapse for both the directed bond percolation model and our own 
using metric factors. Insert: Data collapse for both the directed bond percolation 
model (top right) and our own (bottom left). 



We use the best known DP values of 5 = 0.159464 and u» = 1.733847 pU] and see clear 
data collapses in both models, especially for large t. Using now the obtained values for 
the metric factors, a tt DP = 1-57, oa,dp = 7.51, a tt p = 3137.30 and oa,p = 0.2, we have 
a complete data collapse of all the data from both models as shown in the main plot in 
Figure El indicating that the scaling function for both models is identical. 

Having examined both the critical exponents and the scaling functions, we conclude 
that our model belongs to DP in 1 + 1 dimensions. 

5. First-Order Phase Transitions 

Due to DP having a critical dimension, d c — 4 [18], it is surprising that this model 
belongs to DP in 1 + 1 dimensions only due to the first-order transition in higher 
dimensions. This is not a unique property to this model. A stochastic cellular automaton 
model developed by Bidaux, Boccara and Chate (BBC) |19j is known to belong to DP 
in 1 + 1 dimensions [20J yet also displays a first-order transition in higher dimensions. 

In order to examine the first-order phase transitions, we borrow a technique 
developed by J. Lee and J.M. Kosterlitz [21], which allows one to determine the order of 
a phase transition in an equilibrium system. The method detects a temperature-driven 
first-order phase transition by Monte Carlo simulations in a finite system of volume L d 
with periodic boundary conditions by examining the histogram of the energy distribution 

N(E; (3, L) = AfZ- 1 ^, L)n(E } L) exp(-(3E). (19) 

M is the number of MC sweeps, Z is the partition function and Q is the number of 
states with energy E. For the g-state Potts models with q ordered and one disordered 
state, N has a characteristic double-peak close to T = T c for energy values Ei(L) and 
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Figure 7. Normalised histogram, N'(p) for different population densities in the 
1 + 1 dimensional (left) and 2 + 1 (right) dimensional cases showing the results for 
a continuous and first-order phase transition respectively. 

E 2 (L) corresponding to the ordered and disordered states. The peaks are separated by 
a minimum at E m (L). If we define 

A(E',p,L,N) = -hxN(E',p,L,J\f), (20) 

then at /3 = (3 C {L) defined by A{E^ (3, L) = A(E 2 ; (3, L), 

A(E m , (3, L, AO - A(E\ ; 0, L, Af) = AF (21) 

where AF is the bulk free-energy barrier between the states. At a first-order phase 
transition, AF ~ L d_1 for L ^> £, where £ is the correlation length, whereas AF is 
independent of L at a continuous phase transition. 

Due to our model being out of equilibrium, we continue only in the spirit of the 
above method and examine the histogram for the population density N(p). For a first- 
order phase transition, we expect a double-peaked structure at p and at p + = p_ = 0.5 
due to the corresponding phase coexistence. For p d > p dc we expect N(p ) > N(p + ) due 
to the greater chance of extinction, and likewise for p d < p dc we expect N(p ) < N(p+) 
due to the greater chance of survival. This in fact gives us an excellent method for 
determining the critical point since it will be marked by N(po) and N(p + ) being equal. 
At a continuous phase transition, however, due to the power-law behaviour p ~ t~ s , we 
expect N(p ) = at the critical point. Both expectations are confirmed in Figure [TJ 
where we plot the histogram at the critical points for both the 1 + 1 and 2 + 1 dimensional 
cases. 

As has been mentioned, the phase-coexistence enables the critical value p dc to be 
estimated. Iteratively finding the critical value for large enough t and different values of 
L and then plotting Pd c {L) vs. 1/L in Fig. [HI we are able to obtain, by an extrapolation 
to the intercept with the y axis, the value of the critical point in 2 + 1 dimensions to be 
0.0973 ±0.0001. Unfortunately, no conclusive numerical evidence has yet been obtained 
in 3 + 1 dimensions to give an accurate value for p dc . 
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Figure 8. Plot to obtain an approximation for the critical point in the 2+1 dimensional 
model by extrapolation of the points to the intersection with the y-axis. 



6. Concluding Remarks 

By examining the critical exponents 5 and rj and the scaling function for the survival 
probability, we have shown that the 1 + 1 dimensional version of our model belongs to 
the class of DP. Higher dimensions however do not, surprisingly, belong to this class due 
to the observed first-order phase transitions. In order to examine these transitions we 
used a technique inspired by equilibrium systems to find the value of the critical point 
in 2 + 1 dimensions. 

In order to answer the question about why our model does not belong to DP in all 
dimensions, we must look at how our model differs from others in the class. One dif- 
ference to other models such as the contact process [22] or the Domany-Kinzel [23| [23] 
model is that, here, population growth requires two particles meeting, whereas death 
involves one particle only This requirement results in the existence of the critical pop- 
ulation density, meaning that extinction can occur either because the death rate is too 
large and/or because the population density is too small. The previously mentioned 
BBC model shares this feature since they too observed a critical value of the initial con- 
centration of "alive" sites below which the system decayed to a zero concentration state. 
With this inclusion, it is remarkable that the 1 + 1 dimensional version of our model 
belongs to DP at all and only goes to highlight the robust nature of the universality 
class. 

It seems likely that the unusual behaviour of both our model and that of BBC 
displaying DP behaviour in 1 + 1 dimensions, whilst first-order transitions in higher 
dimensions, is due to the critical population density/concentration present in both 
models. In our model, with the given rules for growth and decay, it is clear that the 
population will become more dependent on the density as the dimensionality of the 
system increases since particles will find it progressively harder to meet one another 
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before they die. Time-dependent simulations carried out in higher dimensions, for 
example, needed increasingly large populations at t — for initial population growth 
as pa was increased, whereas two particles were sufficient in 1 + 1 dimensions even for 
< p d ~Pd c < 1. 

To examine whether the first-order transitions are a result of the critical population 
density, we changed our model to include single-particle reproduction, A — ► 2A, at rate 
c. From mean field analysis, this inclusion eradicates the critical population density and 
initial numerical analysis appears to show that this changes the transition to continuous 
in at least the 2 + 1 dimensional case. 

In the 1 + 1 dimensional case, not only is it easier for particles to meet, but the 
larger fluctuations in the population are likely to be enough to induce the observed 
continuous phase transition as has been known to happen [25]. This was tested by 
introducing two separate types of particles into the system, L and R particles, where 
the L particles could only move left and R particles to the right. All other rules remained 
the same except that both an L and a R particle were required for reproduction. This 
modification to our model reduced the size of the fluctuations in the overall population 
and, again, early numerical simulations appear to show that this results in the transition 
changing to first-order. 

Clearly, more work will however have to be carried out to confirm the above 
observations and to find out exactly what is happening at the phase transition. 

Acknowledgments 

We would like to thank Daniel Lawson and Gunnar Pruessnar for useful discussions and 
Uwe Tauber and Iwan Jensen for helpful email correspondence. Alastair Windus would 
also like to thank the Engineering and Physical Sciences Research Council (EPSRC) for 
his Ph.D. studentship. 

References 

[1] Hayc Hinrichscn. Non-equilibrium critical phenomena and phase transitions into absorbing states. 

Adv. Phys., 49(7):815-958, 2000. 
[2] O. Peters and J. D. Neelin. Critical phenomena in atmospheric precipitation. Nat. Phys., 2(6):393- 

396, 2006. 

[3] A. Vespignani, R. Dickman, M.A. Munoz, and S. Zapperi. Absorbing-state phase transitions in 

fixed-energy sandpiles. Phys. Rev. E, 62(4):4564-4582, 2000. 
[4] S.M. Dammer and H. Hinrichsen. Epidemic spreading with immunization and mutations. Phys. 

Rev. E, 68(1):016114, 2003. 
[5] S. Liibeck. Universal scaling behaviour of non-equilibrium phase transitions. Int. J. Mod. Phys. 

B, 18(3977), 2004. 

[6] P. Grassberger. On phase transitions in Schlogl's second model. Z. Phys. B, 47:365-374, 1982. 
[7] H. K. Janssen. On the non-equilibrium phase-transition in reaction-diffusion systems with an 

absorbing stationary state. Z. Phys. B, 42(2):151-154, 1981. 
[8] I. Jensen and R. Dickman. Noncquilibrium phase-transitions in systems with infinitely many 

absorbing states. Phys. Rev. E, 48(3):1710-1725, 1993. 



Phase transitions in a lattice population model 



12 



[9] H. Hinrichscn. First-order transitions in fluctuating 1+1-dimensional non-equilibrium systems. 
cond-mat/0006212, 2000. 

[10] I. Jensen. Low-density series expansions for directed percolation: I. A new efficient algorithm with 
applications to the square lattice. J. Phys. A-Math. Gen., 32:5233 - 5249, 1999. 

[11] F. Schlogl. Chemical reaction models for non-equilibrium phase transitions. Z. Phys. B, 253:147- 
161, 1972. 

[12] S. Lubeck. Scaling behavior of the order parameter and its conjugated field in an absorbing phase 
transition around the upper critical dimension. Phys. Rev. E, 65(4):046150, 2002. 

[13] P. Grassberger and De La Torre. Reggeon field theory (Schlogl's first model) on a lattice: Monte 
Carlo calculations of critical behaviour. Ann. Phys., 122:373, 1979. 

[14] P. Grassberger. Directed percolation in 2+1 dimensions. J. Phys. A-Math. Gen., 22(17):3673- 
3679, 1989. 

[15] I. Jensen. Critical behaviour of the three-dimensional contact process. Phys. Rev. A, 45:R563, 
1992. 

[16] I. Jensen. Critical exponents for braching annihilating random walks with an even number of 

offspring. Phys. Rev. E, 50(3623), 1994. 
[17] I. Jensen. Low-density series expansions for directed percolation: III. Some two dimensional 

lattices. J. Phys. A-Math. Gen., 37:6899-6915, 2004. 
[18] J.L. Cardy and R.L. Sugar. Directed percolation and Reggeon field theory. J. Phys. A, 13:L423, 

1980. 

[19] R. Bidaux, N. Boccara, and H. Chate. Order of the transition versus space dimension in a family 

of cellular automata. Phys. Rev. A, 39(6):3094-3105, 1989. 
[20] I. Jensen. Universality class of a one-dimensional cellular automaton. Phys. Rev. A, 43(6):3187- 

3189, 1991. 

[21] J. Lee and J.M. Kosterlitz. New numerical method to study phase transitions. Phys. Rev. Lett., 
65(2):137-140, 1990. 

[22] T. E. Harris. Contact interactions on a lattice. Ann. Probab., 2(6):969-988, 1974. 
[23] W. Kinzcl. Phase transitions of cellular automata. Z. Phys. B, 58:229-244, 1985. 
[24] E. Domany and W. Kinzel. Equivalence of cellular automata to Ising models and directed 

percolation. Phys. Rev. Lett., 53:311-314, 1984. 
[25] J. Marro and R. Dickman. Nonequilibrium phase transitions in lattice models. Cambridge 

University Press, 1999. 



